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ABSTRACT 

The effects of stochasticity on the luminosities of stellar populations are an often neglected but 
crucial element for understanding populations in the low mass or low star formation rate regime. 
To address this issue, we present SLUG, a new code to "Stochastically Light Up Galaxies". SLUG 
synthesizes stellar populations using a Monte Carlo technique that treats stochastic sampling properly 
including the effects of clustering, the stellar initial mass function, star formation history, stellar 
evolution, and cluster disruption. This code produces many useful outputs, including i) catalogs of 
star clusters and their properties, such as their stellar initial mass distributions and their photometric 
properties in a variety of filters, ii) two dimensional histograms of color-magnitude diagrams of every 
star in the simulation, iii) and the photometric properties of field stars and the integrated photometry 
of the entire simulated galaxy. After presenting the SLUG algorithm in detail, we validate the code 
through comparisons with SB99 in the well-sampled regime, and with observed photometry of Milky 
Way clusters. Finally, we demonstrate the SLUG's capabilities by presenting outputs in the stochastic 
regime. SLUG is publicly distributed through the website http://sites.google.com/site/runslug/ 
Subject headings: methods: statistical; galaxies: star clusters; galaxies: stellar content; stars: forma- 
tion; methods: numerical; techniques: photometric 



1. INTRODUCTION 

Fundamental progress in understanding the proper- 
ties of galaxies, star clusters and stellar populations 
comes from the comparisons between observed photom- 
etry and synthetic photometry derived from stellar evo- 
lution codes. It has become common practice to infer 
properties such as star formation rate (SFR), star forma- 
tion history (SFH), age, metallicity, redshift, and stel- 
lar mass from photometry. Despite the limits of theo- 
retical modeling of stellar populations (such as uncer- 
tainties with dust, stellar evol ution, and the stellar ini- 
tial mass function (IM F); see iConrov eTaIll2009l [20T0I : 
iConrov fc GunnI 120101 ) synthetic libraries have reached 
a degree of precision that allows accurate estimates of 
these parameters - although sometimes with degeneracy 
- in massive galaxies and clusters. 

However, observations reveal a higher complexity in 
lower mass systems where scaling relations which ap- 
ply to more massive systems cannot be trivially extrap- 
olated (e.g., |LeEet^[2002 [WiEiFin [200^ More- 
over, in lower mass systems, the limited number of stars 
that are present invalidates the basic assumption used by 
most of the currently available codes for syntheti c pho- 
tome try (such as STARBURST99 (SB99: .Leithcrer et al.l 
[1993); P EGASE (iFioc fc Rocca-Volmerangel [l997^: and 
GALEV (jKotulla et al.ll2009| )): that the "IMF is fully 
sampled at all times. Violation of this assumption leads 
to stochastic variations in photometric properties that 
these codes do not fully capture. 

For example in globular clusters, the simplest observed 
stellar populations, failure to account for sampling ef- 
fects can lead to a dramatic overestimate of the contri- 
butions of blue horizontal branch and AGB stars to the 
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integrated light. As a result, correct estimates of glob- 
ular cluster ages and metallicities based on their inte- 
grated light are possible only if one correctly accounts 
for stochasticity |Colucci et al.ll20l"T[ ). 

Moreover, in weakly star forming regions, 
stochastic effects can mimic those of a vary- 
ing IMF. Indeed, recent observations in the low 
SFR regime hav e led to serious consideratio n of a 
varying IMF (IPflamm-Altenburg fc Kroupal 
Hoversten fc Glazebrookl 120081: IMeurer et all 



2008 



2009; 



Lee et al.ll2009[ ). However a fully self-consistent model of 



stochasticity, allowing for a full range of parameters such 
as differing degrees of stellar clustering, metallicities, 
steUar tracks, input IMFs and CMFs, and SFHs has 
not been available to test the null hypothesis of a 
non-varying but stochastically sampled IMF. 

These consid erations ap p ly not only to the dwarf galax- 
ies studied by 'Le e et all (|2009[ ) but also to the outer 
region s of galaxies s uch a s XUV disks (jBoissier et al.l 
I2007al: [Th ilkcr et &\. 20 071) and outlyin g H II regions 



( Werk et ah 2008.: Gogarten et aLll20Q9D where stochas- 
ticity becomes crucial in the interpretation of inferred 
SFRs and SFHs. 

While the number of studies that use Monte Carlo ap- 
proaches to address problems on scales of clusters and 
galaxies is growing, a general purpose tool to study pho- 
tometry in clusters and galaxies has not previously been 
available. To fill this need, we have created SLUG, a 
code to allow proper study of the stochastic star forma- 
tion regime at a range of scales from individual star clus- 
ters to entire galaxies. SLUG provides a variety tools 
for studying the stochastic regime, such as the ability 
to create catalogs of clusters including their individual 
IMFs and photometric properties, color-magnitude dia- 
grams (CMDs) of entire galaxies where we keep track of 
the photometry of every star, as well as integrated pho- 
tometry of entire composite populations. 
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This paper, the first of a series, focuses on the raethods 
used in the code along with several tests to demonstrate 
that we are reliably reproducing observations and other 
synthetic photometry predictions. We then demonstrate 
the use of this cod e in the stoc hastic regime. In a com- 
panion paper (Fumagain et all [20 Hi ), we use SLUG to 
demonstrate that, once random sampling is included, a 
stochastic non-varying IMF can reproduce the observed 
variation of the Ha/FUV ratio in dwarf galaxies without 
resorting to modifications of the IMF. In a the second pa- 
per of the series (da Silva et al in prep.) we will explore in 
detail the implications of stochastic star formation with 
clustering. Further work will apply this code to a variety 
of astrophysical questions, such as understanding SFR 
calibrations in the stochastic regime and further study 
of other claims of a varying IMF. 

The layout of the paper is as follows: 321 presents an 
introduction to stochasticity and its effects on the lumi- 
nosity of stellar populations; SJS] gives a detailed descrip- 
tion of the SLUG algorithm; gl discusses various tests of 
the code; !j5]shows a presentation of the code's outputs in 
the stochastic regime; finally, !j6] summarizes the results. 

2. WHAT IS STOCHASTICITY? 

Many astrophysical studies require creation of syn- 
thetic photometry of galaxies and other collections of 
stars in order to compare with observations. In this 
section we present a discussion of the various effects of 
stochasticity and the regimes in which they are impor- 
tant. 

2.1. Coeval Stellar Populations 

The standard procedure for calculating the luminosity 
from a coeval population of stars used by the most pop- 
ular implementations (such as SB99) is as follows. To 
find the luminosity per unit mass of a coeval population 
in some band /3 at a time t after formation (i/s^coevaiit)), 
one simply integrates the luminosity per unit mass of 
each star in that band as a function of mass and time 
{lp{m,t)) weighted by the distribution of stellar masses 
(i. e. the IMF) 

^l3,coevai{t) ^ I ^f}{'m,t)— dm. (1) 

Jm,ni„ a mm 

Note that here we use a normalization of the IMF such 
that f"""^ TTT^dm = 1. 

^^min alnm 

By performing this integral, these models assume an 
infinitely well-sampled IMF. As a result the above for- 
mula is mass-independent, meaning that ({S^coevai can be 
scaled according to the total amount of stellar mass in 
a population (i.e. the luminosity of a mass M of stars 
is simply Mi/j^coevai)- Thus a given amount of mass M 
will have a 1-to-l mapping to a particular luminosity L. 
However for small stellar populations, the assumption of 
continuous sampling breaks down and effects of stochas- 
ticity can become important. Specifically, stochastic ef- 
fects create a statistical dispersion of luminosities that 
result from a given mass M of stars based entirely on the 
probabilistic sampling of the mass distribution of stars. 
This is because each realization of a given mass M is built 
up with a different distribution of stellar masses which, 
due to the non-linear dependence of luminosity on stellar 



mass, yields a different luminosity. We call this type of 
stochastic process sampling stochasticity. 

Perhaps the most important manifestation of sampling 
stochasticity is undersampling of the upper end of the 
IMF. Since the IMF is steeply declining with increasing 
stellar mass, the expectation value of a low mass popu- 
lation drawing a massive star is small. As a result, the 
IMF in a low mass population with few stars can appear 
truncated and have less luminosity than a fully-sampled 
assumption would have predicted. This is due to the very 
super linear dependence of luminosity on stellar mass. 

One can roughly estimate the mass below which this 
effect is insignificant by calculating the expectation value 
of obtaining a star above a given mas s. We do so follow- 
ing the formalism of lElmegreenI (|2000f) , who find that the 
total mass {M) required to expect a single star above a 
mass m is 

/ X 1.35 

Af ^ 3 X 10^ ( ] . (2) 

\100MqJ ^ ' 

This statement i s clear ly dependent on one's choice of 
IMF. lElmegreenI (|2000D uses a Salpeter IMF with a a 
lower limit of 0.3 Mq and no upper limit. If one im- 
poses an upper limit to the stellar mass function, this 
relation turns over and asymptotically approaches the 
limit. However, for order-of-magnitude purposes here, 
we neglect such consideration. This result implies that 
in order to reasonably expect even a single 120 Mq stai0, 
one would need at least a total mass sampled of approx- 
imately IO^'M© = Mtrunc- Thus this IMF truncation 
effect of sampling stochasticity can be ignored for co- 
eval populations with masses 3> Mtmnc- For more ref- 
erenc e on the limits of stochastic sa mpli ng, we recom- 
mendJCervino fc Vails- GabaudI ()2003D and lCervifio et al.l 
(j2003[) . For specific considerations to Ha luminosity (one 
of the features of a stellar population most s ensitive to 
stochasticity), see lCervino fc Luridianal ([2004') . 

Another manner in which stochastic sampling can 
manifest in coeval populations is for stars going through 
particularly short-lived and luminous phases of evolution 
after they leave the main sequence (e.g., AGB and blue 
horizontal branch stars). Since these phases are short, 
only a very narrow range of masses is undergoing one 
of them at any given time. Thus the exact sampling 
within that mass range can have a large impact on the 
number of stars within that phase. As a result, a non- 
infinite population of stars can have additional random 
scatter in luminosity even if M > Mtmnc- This effect is 
more important in populations with little ongoing star 
formation relative to their stellar mass (otherwise new 
stars dominate the photometric properties of the popu- 
lation), at specific ages when these post- main sequence 
populations con tribute significantly to the luminosity of 
the population (jColucci et al.l[20Tll ). 

2.2. Composite Stellar Populations 

In order to characterize a more complicated star for- 
mation history, SB99 and other such schemes integrate 

Due to limitations of stellar evolutionary tracks, this is the 
highest stellar mass SLUG can model and is a reasonable guess 
for the highly uncertain absolute stellar mass limit. While some 
e.g.. [Pi ger 2005) suggest a value of ~ 120 — ISOMq others 
Crowther et al. 2010i 'l suggest it may be as high as 300 Mq. 
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Fig. 1. — A schematic flow-chart describing the algorithm of the the SLUG code. Note that for the case of unclustered star formation, 
the cluster mass is drawn from the IMF and the population step is skipped as the single star is treated as part of a disrupted cluster for 
the remainder of the code. Note this is updated from Fumagalli et al. (2010'). 



over the coeval populations discussed above to find the 
luminosity of all stars in a given band at a time r, 

LfS,total{T) = SFR{t)ei3^coeval{T - t)dt. (3) 

J — OQ 

Such a treatment makes two key assumptions: (1) each 
of the summed coeval populations is large enough to ig- 
nore the effects of sampling stochasticity and (2) the SFR 
is continuously sampled as well. These assumptions can 
quickly break down for sufficiently low SFRs. 

To illustrate this point, consider a galaxy forming stars 
at a constant rate. In order for the IMF not to be trun- 
cated within some time interval dt, there need to be at 
least Mtrunc worth of stars formed in that interval. For 
the SFR to be considered reasonably well sampled, dt 
must be much smaller than the evolutionary timescales 
of any of the stars, which are « 10^ yr for the massive 
stars that generally dominate the light in an actively star- 
forming system. Thus these assumptions require 



SFR 



< lOVr. 



(4) 



Thus these effects can only be ignored for SFRs consis- 
tently » lO^^A/0yr^^ = SFRjemp- However, this tem- 
poral stochasticity is amplified when one considers that 
stars are believed to be formed in discrete collections 
known as clusters. As a result, the clumping in time of 
star formation in clusters can produce stochastic effects 



even in regions with SFRs higher than SFR, 



temp ■ 



In this 



case the characteristic mass in Equation |3] is replaced 
with a mass characteristic of the clust ers being drawn 
(discussed further in da Silva in prep.; iFumagalli et al.l 
[2011.) . 



The conditions required to treat a stellar popula- 
tion as continuous (as opposed to stochastically sam- 
pled) break down in a variety of a strophysical env iron- 
ments such as dwarf galaxies (e.g. jLee et al.|[2009t) , low 
star formation rate regions in the outskirts of galaxies 
(e.g., [B oissicr eTaD l2007bl : IFumagalli fc Gavazzil [200l 
iBigiel et al.ii201 0^. and low surface brightness galaxies 
(e.g.. lBoissieret al.. .200&') . 

3. TECHNIQUE 
3.1. Overview 

Here we present a brief overview of the code while we 
present each step in detail in the subsequent sections. 

SLUG simulates star formation according to the 
scheme presented in Figure [TJ We create collections of 
star clusters obeying a user-defined cluster mass function 
(CMF) (which can include a given mass fraction of stars 
not formed in clusters), SFH, IMF, and choice of stellar 
evolutionary tracks, which we call a "galaxy" . A descrip- 
tion of the parameters that users can vary is provided in 
Table [H 

These galaxies are built up ( §3.2p by first drawing the 
mass of an individual cluster from a CMF. This cluster's 
mass is then filled up with stars according to an IMF. The 
age of the cluster is drawn from a distribution weighted 
by the given SFH. Each of the stars within the cluster is 
evolved using a stellar evolutionary track combined with 
a model spectral energy distribution (SED) to determine 
a variety of integrated fluxes corresponding to commonly 
used photometric filters ( H3.3p . 

At a given set of time steps, these fluxes are summed 
over each star cluster. The clusters are then disrupted 
according to the prescription of .Fall et al- (,2009a ). Dis- 
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TABLE 1 
Input Parameters 



Parameter 


Description 




Controlling the Physics 


IMF 


stellar initial mass function; can choose 




Kroupa, Salpeter, Chabrier, IGIMF, or 




an arbitrary slope 


CMF 


cluster mass function, can change 




slope, minimum and maximum mass 


Stellar Evolutionary 


library of models used for stellar evolution 


Tracks 




Metallicity 


metallicity of the stellar population 


Stellar Atmosphere 


which scheme and models are used for SEDs 


Stellar Wind Model 


which wind model is used for SEDs 


Fraction of stars in 


mass fraction of stars formed in clusters 


clusters 




Controlling the Simulation 


Maximum time 


how long the simulation is run 


SFH 


can be arbitrary 


Seed 


random seed used for simulation 


Controlling Output 


Time step 


time between code outputs 


Fluxes 


choose which fluxes to output 


Colors 


which colors to use for CMDs 


CMD output 


choice of number of bins and 


parameters 


range of color and luminosity for each CMD 


Cluster output? 


set to print output for each cluster 


IMF output? 


set to output IMF histograms for each cluster 




Fig. 2. — IDL GUI interface for running the code. The code 
may also be called via the UNIX or IDL command lines. 



rupted clusters have their fluxes added to a "field" pop- 
ulation while surviving clusters have their properties 
stored individually. The code repeats this process un- 
til a stellar mass equal to the integral of the provided 
SFH is created. 

The code outputs a variety of files that keep track of 
the properties of the stars, clusters, and total integrated 
stellar populations. Table [2] provides a short description 
of each available output file. All outputs are parsed and 
transformed into binary FITS tables. 

The code is open source and written in C-I--I- with 
wrapping and parsing routines written in IDL. This en- 
tire process can be controlled through an IDL graphi- 
cal user interface (see Figure [2]) or either of the UNIX 
or IDL command lines. The IDL routines are wrapped 
in packages for use with the IDL virtual machinal 
for those without IDL licenses. For a full manual 
on how to use the code, visit the SLUG website at 
http: / /sites. googl e . com/site/runslug / 

3.2. Cluster Creation 

Most stars are th ought to be born in star clusters 
(|Lada fc Ladaf |2003l ) and the distribution of star cluster 
masses appears to o bey a power law dist ri bution, where 
observations (e.g.. iZhang fc Falll 119991: ILada fc Ladal 
[2nnl iFall et al.ll2n09hl: IChandar et al. [|20100 and theory 




20 40 60 80 100 

Age [Myr] 



Fig. 3. — Examples of star formation histories average over 1 Myr 
bins for simulations with varying input constant SFRs of 0.0001- 
100 Mq yr~^. The dotted lines show the input SFR. The average 
SFR of the simulation in each case is within 2, 0.2, and <0.02 per- 
cent of the input for 10~*,10~^, and > 10~^ Mq yr~^ respectively. 
SFRs of zero are masked. 



^ which is available for free from 

|http://www.ittvis. com/language/en- us7productsservices/idl/idImodules/idlvirtualmachine.aspx| 
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TABLE 2 
SLUG Output Files 



Name 



Description 



Histogram a 2d histogram of the user's choice of color-magnitude diagram(s) 

of every star in the "galaxy" at each timestep 
Cluster mass, fluxes, most massive star, number of stars, and 

age of each undisrupted cluster at each timestep 
IMF a histogram of the IMF of each cluster that appears in the Cluster file 

Integral the total flux of the entire "galaxy" at each timestep 

Miscellaneous the total stellar mass actually formed, 

as well as the actual SFH and CMF of the simulation 



(e.g-. lFall et~aIll2010D suggest that the index (/3) of the 
power law dN/dM oc is approximately 2. SLUG 

allows for both clustered and unclustered star formation. 
The user can choose what fraction of all stellar mass they 
wish to form in star clusters. If the code is forming clus- 
ters, the CMF's power law slope as well as its upper and 
lower bounds can be varied. If unclustered star forma- 
tion is desired, the stars' masses are drawn individually 
from an IMF and treated as a disrupted "cluster" of one 
star for the remainder of the code. 

The initial masses of stars are drawn from a, n 
IMF. C h oices of IMlgl c u rrently are iChabrieii (|2003D . 
iKroupal ()200in . iSalpeteii (jlQSSj) . a user-defined arbi- 
trary power law, and the recently proposed IGIMF 
(IKroupa k, Weidne r 2003; P flamm-Altenburg fc Krotipal 
120081). While the Chabrier. Kroupa, Salpeter, and power 
law IMFs are implemented as a standard probability den- 
sity function of stellar masses, the IGIMF has additional 
features that require different treatment (see Appendix 

Regardless of the choice of IMF, we draw stars until the 
total mass of the star cluster is built up. Since the ran- 
dom distribution of stars never exactly equals the mass 
of the cluster, a question arises as to whether to keep the 
last star added. This last star increases the mass of the 
cluster above the cluster mass drawn from the CMF. We 
determine whether or not to keep that star in the clus- 
ter based on whether keeping the star in makes the total 
mass of stars closer to the mass drawn from the CMF 
than leaving it ou10. 

Independent of its mass, the age of the cluster rela- 
tive to the galaxy is assigned in a probabilistic manner 
weighted by the SFH (which can be arbitrary) such that 
the SFH is reproduced on average. This produces a scat- 
ter in the SFHs for even a given "constant" SFR. Thus 
slug's definition of a galaxy with a constant SFR is 
not a galaxy where the SFR is constant at every individ- 
ual tim^U, but rather a galaxy that produces an amount 
of stars over a time dt equal to SFRxdi which are dis- 
tributed in clusters whose ages are drawn from a uniform 
distribution. This interpretation of what a SFR is and 
its implications is discussed in more detail in da Silva et 

^ IMFs are truncated at 0.08-Mq due to lack of lower mass stellar 
tracks 

The effects of different sampling methods and thei r depen dence 
on the CMF is studied in detail bv IHaas fc Anders! l|20ia ). Our 
method is identical to their 'stop- nearest' method. 

^ A constant SFR cannot be instantaneously constant because 
stars form in discrete units of mass. For example, when a star is 
born, the instantaneous SFR is infinite, thus we must turn to a 
more probabilistic interpretation of the SFR. 



al. (in prep.). 

Clusters are born until the total mass of stars formed 
is equal to the integral of the SFH. As with the problem 
of populating a cluster with stars, a galaxy will never be 
filled to exactly its given mass with an integer number 
of clusters. Therefore we apply the same condition for 
populating the galaxy as we do the clusters: we add until 
we exceed the mass and keep the final cluster only if the 
total galaxy mass is closer to the desired value if we keep 
it. As a result the average SFR over the entire simula- 
tion of a particular galaxy can be higher or lower than 
the input value. This effect is small for most regimes, 
but very rare drawings of the CMF at low SFRs can pro- 
duce mild departures. We emphasize that this is not the 
effect of any error associated with the code but rather is 
the necessary result of our interpretation of what a SFR 
means. 

We demonstrate the results of this procedure in Figure 
[31 The figure shows that, while lower average SFRs tend 
to produce larger fractional scatter in the instantaneous 
SFR, significant scatter remains until SFRs exceed 10 
Mq yr^^. This scatter is a direct result of the finite size 
of clusters. To clarify with an example, consider that a 
1O''M0 cluster (when averaged over the 1 Myr similarly 
to the curves shown in Figure [3]), will appear as a deviant 
peak for all but the highest SFRs, where the contribution 
of that individual cluster is drowned out by enough other 
clusters. 

We note that in this release of the code all stars in a 
cluster are treated as having identically the same age. 
While observations suggest a scatter o f a several Myr 
(IPalla fc Stahlei< [19991 i Jeffriel [2007t IHosokawa etall 
[2011D, the mass dependence of this scatter is unclear. 
Given the uncertainties, and that the intracluster age 
scatter is at most a few Myr, we chose to neglect this 
effect for now but plan on implementing it in the future. 

3.3. Stellar Tracks, SEDs, and Broad Band Photometry 

Given the mass and age of each star, we need to de- 
termine its properties for a variety of observables. Our 
method uses many of the same algorithms found i n SB99 
(jLeitherer et al.ll 19991: IVazquez fc Leithereiil2005D to cre- 
ate a set of tables from which SLUG interpolates. These 
tables are constructed in advance so they need not be 
computed at run time. 

Our first step is to determine the physical properties of 
each star. To this end, we make use of a variety of stellar 
evolutionary models. Modifying the SB99 source code, 
we were able to obtain the full range of stellar tracks 
available to SB99 (see Table [3]). In the future we plan 
to implement a wider range of stellar tracks including 
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those f romlEldridge fc Stanwa"\^ (I2009D and th e BaSTI h- 
brary ()Pietrinferni et aLl 120071: ICordier et aLll2007[ ). We 

supplement the Geneva tracks with the Padova+AGB 
tracks for stars in the mass range 0.15-0.8 Mq. These 
models provide luminosities, gravities, chemical compo- 
sitions, and effective temperatures at discrete intervals 
in the evolution of a discrete number of stellar masses. 
We then need to map these physical properties to stellar 
atmospheres in order to estimate the spectral energy dis- 
tributions of the stars. Our code allows users to choose 
from one of five possible SB99 algorithms for modeling 
the atmospheres. We implement all four prescriptions of 
stellar winds available in SB99 (Maeder, empirical, theo- 
retical, and Elson), which affect the SEDs for Wolf-Rayet 
stars for some regimes and prescriptions. It is important 
to note that the SB99 algorithms match SEDs to tracks 
with a nearest neighbor approach and not through in- 
terpolation. Therefore there can be some discreteness in 
the output SEDs. Future work will include removal of 
this effect. 

With SEDs in hand, we can convolve with filters to 
determine the photometry of each point in our stellar 
tracks. For this step we include the effects of nebular con- 
tinuum (free-free, free-bound, and 2 photon processes) as 
implemented in SB99, but neglect nebular line emission 
for this first release of the code. (For a discussion of 
the importan ce of nebular c ontinuum for the SEDs, see 
iReines et al. , 2 10[ A lso see iLeitherer fc HeckmanI 119951 
and lMoUa et al ."2009'.) The fuU list of available fihers is 
presented in Table ID We also integrate the SED to de- 
termine the bolometric luminosity as well as to calculate 
Q(H''), the number of hydrogen ionizing photons emit- 
ted per second. One can convert Q(H*') to Ha luminosity 
with a simple convers ion assuming case B re combination 
(our notation follows lOsterbrock fc Ferla"iidl[2C)06il . 

Lhc = (1 - fesc){l - /d«.t)Q(H°) huHc. 

« 1.37 X 10-12(1 _ /^^^)(i _ /rf„,,)Q(H°) ergs/s 

(5) 

where fesc is the escape frac tion (tho ught to be be- 
tween 0.05 (jBoselh et al.ll2009[ ) and 0.4 (jHirashita et al.l 
l2003f )) and fdust represents the fraction of of ionizing 
photons ab sorbed by du st grains (e.g., see appendix of 
iMcKee fc W ilhams ll997L who suggest a value of 0.37). 
To better characterize the ionizing luminosity we also 
keep track of Q(He'') and Q(Hei) which represent the 
numbers of ionizing photons in the He I and He II con- 
tinua respectively. 

The above steps allow us to create a discrete two- 
dimensional table for each flux band where one axis rep- 
resents stellar mass, the other represents time, and the 
value of the table is the logarithm of the fiux in that band 
at the appropriate mass and time. Our tables are created 
through use of the isochrone synthesis method such that 
our results are stable against the numerical issu es that 
arise from a fixed mass approach (Chariot fc Bruzuall 

[MiD. 

3.4. Evaluating the Stellar Properties 

To determine the properties of a given star of any mass 
at any given time, we first determine if the star is still 



alive. This is done by an interpolation in time to find 
the minimum mass of a dead star {rrideath) at a given 
time according to our stellar evolution models (where 
we call a star "dead" if it no longer has entries in our 
stellar tracks). If the star is less massive than nideath, we 
interpolate our model tables to determine the fiux in a 
given filter within 0.01 dex . 

For computational speed, there are a variety of approx- 
imations and restrictions we are forced to implement. 
The current scheme only allows ages up to 1 Gyr for 
the stellar tracks (to be expanded in later releases of the 
code). We do not evolve stars less massive than 0.9 Mq 
(a number which can be changed by the user). These 
stars do not evolve past the main sequence for the cur- 
rent maximum age of the code of 1 Gyr, so these stars are 
treated as having their zero-age main sequence (ZAMS) 
properties. Due to limitations of the stellar tracks, we 
treat the photometric properties of all stars less massive 
than 0.156 Mq identically to those of 0.156 Mq stars. 
For many purposes, more massive stars dominate the 
light in the bands such that this approximation is rea- 
sonable. The tracks also impose a 12OAf0 upper mass 
limit on stars. 

Curre ntly, we neglect the eff ects of binary stellar evolu- 
tion (see lEldridge fc Stanwav.2009 ). which may have an 
impact on the derived results by producing a bluer pop- 
ulation with a reduced number of red supergiants and 
increased age range of Wolf-Rayet stars. 

3.5. Cluster Disruption 

If the user chooses to form stars in star clusters, we 
randomly disrupt our clusters in a mass independent way 
such that dN/dr cx r'^ (following .Fall et al. 2009a). We 
start cluster disruption 1 Myr after the cluster forms. 
This results in 90% of star clusters being disrupted for 
each factor of 10 in age after 1 Myr. Stars in disrupted 
clusters still have their photometry calculated for the in- 
tegrated properties of the galaxy and are kept track in a 
set of "field" variables and outputs. 

4. VALIDATING TESTS 

In this section we present a variety of tests to validate 
the outputs of SLUG. For these tests we make use of a set 
of fiducial pa rameters presented in Table [5] unless other- 
wise notecjrl. To emphasize that SLUG can be applied 
at different scales, we arrange these tests in order of scale 
starting with individual clusters and then considering in- 
tegrated properties of entire galaxies in the well-sampled 
regime. 

4.1. Photometry of Clusters 

To demonstrate that SLUG reproduces properties of 
observed clusters, we turn t o the c atalog of young star 
clusters compiled in ILarsenl (jl999() . To reproduce the 
clusters we modify our fiducial IMF to extend down to 
0.08 Mq and run a SLUG model with a SFR of IMq 

* While the preferred SEDs for SB99 are the Pau-(-Smi atmo- 
spheres, we find that the Pauldrach models are far too discrete. 
Therefore while we provide the Pau-f Smi atmospheres, we recom- 
mend the Lej-(-Smi. 

^ Since we aim to test SLUG rather than to perform a study 
of the effects that the multiple parameters have on the luminosity 
distributions , we choose widely adopted vlaues. 
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TABLE 3 




Stellar, Properties 


Parameter 


Allowed Values 


Tracks 


Geneva STD^, Geneva High"-, Padova STD*', Padova AGB'' 


Metallicity<= 


0.0004-0.50 


SEDs 


Planck'^, Lejeune'=, Lejeune+Sch^ , Lejeune+SMF, Pau+SMl'^ 


Wind Models 


Maeder', Empirical', Theoretical', Elson' 



"[Mcynct ct al. (1994) and references therein 
^ [Fagotto et ah, (.1994.) and references therein 
" solar is 0.20 

simple blackbody SE P 
^ ILeieune et all 09% [1998) 

^ same as e, but for stars with strong winds use lSchmutd 11998) 

s same as e, but for stars with strong winds use IHiUier fc Miiiea SJMB ) 

same as g, but use IPauldrach et al.l l l200ll ) for O stars 
' ILeitherer et al.1 STM^ 



> 



pq 

I 




-6 



10' 
Age [Myr] 



0.0 0.5 
V-I 



1.0 



1.5 
1.0 

0.5 

0.0 
-0.5 




"•• ■::=■• :'::-i'v;;:'/;f!;.'*.'i.\ 



> 



0.4 
0.2 
0.0 
-0.2 
-0.4 



■'■'.'■MrV-- 




■./.-i»,!-'Va'.-.-: ■ 



-0.4 -0.2 0.0 0.2 0.4 
B-V 



-0.5 



0.0 0.5 
V-I 



1.0 



Fig. 4. — Comparison of observed young star clusters from [LarsenI (119991 ) (black points) to SLUG models of clusters > 10^ Mq (blue 
triangles). T he or ange curve shows the trajectory of a SB99 10^ cluster. Data are omitted from upper left panel as the ages are not 
prese nt in the lLars en ( 1999) catalog. Arrows denote the extinction vector for Ay = 0.5 mag (created following appendix B of .Schlegel et al.l 



yr~^ for 500 Myr, evaluated every 10 Myr. Note that 
the the SFR does not directly affect the CMF or the 
properties of the clusters, only the number of clusters in 
existence at a given time. We show the results of this 
exercise in Figure |4] where we find remarkable agreement 
between the models and the data. As is clear from the 
figure, we are able to reproduce both the location and 
spread of most of the observed data. Clusters that fall 
outside of the locus of SLUG models fall can easily be 
reproduced when one accounts for a modest amount of 
reddening (see reddening vector). 



4.2. Cluster Birthline 

Another test of the photometry of clusters is to com- 
pare their H a luminosity to the ir bolometric luminos- 
ity. Work bv lCorbelli et al.l (|2009( ) has shown that newly 
born clusters lie along a birthline in this parameter space. 
In Fig. [5] we compare the same models as Section |4TT] (as- 
suming fesc and fdust = 0) with their observational data 
and find good agreement. Our theoretical predictions 
di ffer slightly in the tilt of the locus of points from those 
bv lCorbelh eFall (l200l . since we characterize the prop- 
erties of our stars in a different manner (making use of 
stellar tracks rather than fitting formulae). To better 
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Fig. 5. — (iop) Here we present the birthline as first discussed by ICorbelli et all II2009I ). Black points and crosses are data from that 
paper with circle-crosses denoting their 'clean' sample. Blue data points are clusters from SLUG. We see that our models are in relatively 
good agreement with observations, {bottom) We present overlays demonstrating the average IMFs in each region of the birthline plot. 



demonstrate the origin of the birthhne we also make use 
of slug's abihty to keep track of the IMF of each indi- 
vidual cluster (see bottom panel of Figure [5]). Here we 
can see that the birthline from left to right forms a se- 
quence of progressively more well-sampled upper ends of 
the IMF. Extremely rare deviants exist below the birth- 
line where more extremely massive (> lOOA/©) stars are 
drawn than average, resulting in being born below the 
birthline. Note that these rare clusters consisting of 



essentially isolated O stars have al so be en reported in 
the Milky Way (de Wit et al. 20M [200l and the SMC 
(jOev et al.l l20G4t iLamb et al.i i2010[ ) in numbers consis- 
tent with stochastic sampling of the IMF. 

4.3. Comparison with SB99 

A third obvious comparison for SLUG is SB99 itself. 
Since SB99 is widely used, it serves as a benchmark for 
SLUG. Indeed, one of the motivations for making use of 
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TABLE 4 
Broad Band Filters 



Filter 


Reference 


NUV 


1 


FUV 


1 


u 


2 


9 


2 




2 


i 


2 


z 


2 


J 


3 


H 


3 


K 


3 


U 


4 


B 


4 


V 


4 


R 


4 


I 


4 


Q(H0) 


5 


Q(HeO) 


5 


Q(Hel) 


5 


^Bol 


6 



1 Morrisscv ct al. ( 200^) 

2 Fukugita et al. ( 1991) 
^ Skrutskie ct al. (20(3) 

* ADDcnzcUcr ct aL. (ji"99l ) 
^ Obtained by integrating SED 
blueward of 912, 504, and 208 
A for QCH"), Q(HeO), Q(Hei) 
respectively. 

^ Given by stellar evolutionary 
tracks. 



TABLE 5 
Fiducial Inputs 



Parameter 



Fiducial Value 



Time step 
Maximum time 
IMF 
CMF 

Stellar Evolutionary Tracks 
Metallicity 
Stellar Atmosphere 
Stellar Wind Model 

Fraction of stars in clusters 



lO'^ yr 
10^ yr 
1-120Mq; slope=-2.35 
20 - WMq; slope=-2 
Padova+AGB 
Solar; Z = 0.20 
Lej+Smi 
Maeder 
100% 



the SB99 tracks and SED matching algorithms is that 
our code should be able to exactly reproduce SB99 if we 
select input parameters that place us in the continuously- 
sampled regime. To that end we now present a variety of 
tests where we compare to SB99 to demonstrate that we 
can reproduce their results in the this regime (the regime 
of a large galaxy-sized amount of stars). 

To compare the outputs of both SB99 and SLUG, 
we choose an instantaneous burst of star formation to 
demonstrate the matching of the codes in both ampli- 
tude and time. We run a SB99 model similar to our 
fiducial model (i.e. IMF slope of -2.35 from 1-120 M©, 
solar metaUicity, Padova-|-AGB tracks, Lej-|-SMI SEDs, 
and Maeder stellar wind models). To meaningfully com- 
pare with SB99 we must choose SLUG input parameters 
such that we are evaluating a population where SB99's 
approximations are valid. We therefore draw a very large 
instantaneous population of IO^A/q. To nullify any pos- 
sible effects of our procedure for populating the clusters, 




20 30 
Age [Myr] 

Fig. 6. — A comparison of SLUG and SB99 simulations of an 
instantaneous burst of IO'^^Mq. We find good agreement between 
the two in both the absolute normalization of the fluxes as well as 
the time-dependent behavior. FUV and r band fluxes are presented 
in units of ergs/s/Hz while Q(H'') is in units of photons/s. In the 
panel, one can see the effects of the discrete SED matching 
techniques implemented by SB99 in the age ranges of 12-18 Myr 
and 20-22 Myr. 



i ^ 



50 Myr 



200 Myr 







500 Myr 



— -X - 
800 Myr 



5.0'10' l.O'lO* 1.5'10'' 2.0'10* 2.5'10'' 

Wavelength [A] 

Fig. 7. — A comparison of SLUG and SB99 photometry for an 
instantaneous burst of 10^ Mq, evaluated at the ages indicated in 
each panel. The grey solid line represents the output spectrum 
from SB99 for such a population. The filled color circles show 
the SB99 integrated fluxes for the filters available to SLUG. The 
black X 's mark the SLUG photometry for the well-sampled model 
described in section [4.31 

we ensure all clusters are very large by modifying the 
fiducial CMF to a restricted range (10*^ - 2 x IO^Mq). 
We present the results in Fig. [6l It is evident that we are 
accurately able to reproduce SB99 in the well-sampled 
regime for integrated "galaxy" properties. We match 
both the amplitude and time evolution in all photometric 
bands. 

This can also be seen by looking at the full SEDs. In 
Figure [3 we present photometry for all 15 of the flux 
bands available for SLUG and compare with the spectra 
and integrated photometry produced by SB99 at a vari- 
ety of time steps. Again we are able to fully reproduce 
the photometric properties in the well-sampled regime 
from FUV to K-hand. 

In both these demonstrations, SLUG matches SB99 
within 0.026 dex for all fluxes at all times. 

5. STOCHASTICITY IN ACTION 

Having demonstrated that SLUG can reproduce real- 
istic clusters as well as reproduce SB99's results, we now 
present outputs of SLUG in the stochastic regime. 
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Fig. 8. — CMDs of an instantaneous 10^ A/© burst population at 
the ages indicated in each paneL Only stars more massive than 
0.9 Mq are binned in the CMDs. The dotted lines show the cor- 
responding theoretical isochrones. The SLUG CMD has been con- 
volved with circular top-hat PSF solely to improve visibility. The 
color bar denotes the number of stars in that region of the diagram. 

5.1. Effects on Coeval Populations 

Recent studies fe.g., lDalcanton et al.ll2009() have shown 
the wealth of information that can be obtained using re- 
solved CMDs of stars within a galaxy. For comparison 
with such studies in the stochastic regime, SLUG pro- 
duces binned 2 dimensional histograms that keep track 
of the user's choice of color magnitude diagrams. Such 
diagrams allow us to directly characterize the effects of 
stochasticity in a coeval population. In Figure [8l we 
compare CMDs produced by SLUG for a lO^M© instan- 
taneous to 

the theoretical isochrones from which they are pro- 
duced. Aside from demonstrating we accurately repro- 
duce the tracks, we are able to see the effects of stochas- 
ticity in leaving rapid phases of evolution unpopulated. 
Note that SLUG is capable of producing such diagrams 
for any given SFH. 

5.2. Effects on Composite Populations 

While individual clusters of stars can be treated as co- 
eval, larger systems are intrinsically built of composite 
populations. One of the most basic composite popula- 
tions one can consider is a galaxy forming stars at a con- 
stant star formation rate. As discussed in Section 12. 2[ 
the value of the SFR will have a significant impact on 
the effects of stochasticity. 

To demonstrate the differences that stochasticity 
makes, we compare SLUG realizations to those of a well- 
sampled SB99 model. In Figure [HI we first examine the 
luminosities for SFRs of 1, 10'^, and lO'^ Mq yr'^ 
with our fiducial values for the CMF and cluster mass 
fraction. For each SFR, we show the mean and median 
of the SLUG runs along with the 5 and 95 percentiles. 

One can clearly see an increase in fractional scatter as 
one decreases the SFR, which can be attributed to the 
more bursty SFHs which are a result of the grouping of 
age in massive clusters. This scatter appears at higher 
SFRs than predicted by our naive discussion in Section 
I2.2l as a direct result of the clustering. In fact, nearly all 
of the scatter seen in Figure [S] is a result of the clustering 
rather than sampling of individual clusters. This is most 
clearly demonstrated by Figure [10] which shows similar 
simulations but with completely unclustered star forma- 



tion. Without clustering the 10~^ Mq yr~^ models have 
approximately an order of magnitude less scatter in the 
5-95 percentile range of the log of the luminosity. We 
see that the unclustered stochastic effects behave as pre- 
dicted in Section [521 where the fractional scatter is small 
, for SFRs ^ 10"^ Mq yr~^ and quickly increases as th e 
SFR decreases (also discussed in lFumagalli et alll2011[ ). 

For a demonstration of the effects of clustering, we 
present the tracks of a subset of individual stochastic 
realizations of clustered star formation in Figure [TT] One 
can see that the Q(H'^) curves are less uniform than the 
R luminosity. This is a direct result of the sensitivity 
of Q(II°) to the youngest, most massive stars. One can 
also see that the scatter increases with decreasing SFR 
as expected. This is to be further discussed in da Silva 
et al. (in prep.) where we elaborate on the effects of 
stochastic star formation when one includes clusters. 

6. SUMMARY 

We introduce SLUG, a new code that correctly ac- 
counts for the effects of stochasticity (with caveats dis- 
cussed in the text) by populating galaxies with stars and 
clusters of stars and then following their evolution using 
stellar evolutionary tracks. Cluster disruption is taken 
into account and a variety of outputs are created. 

We present a series of tests comparing SLUG to obser- 
vations and other theoretical predictions. SLUG is able 
to repr oduc e the photometric properties o f clusters from 
the jLarsenl (i999f ) catalog as well as the ICorbelli et al.l 
lj2009f ) birthline. It can also reproduce the results of SB99 
in the well-sampled regime. 

Finally we present SLUG outputs in the stochastic 
regime and demonstrate the flexibility of the code to ad- 
dress a variety of astrophysical problems with its variety 
of possible outputs. 

SLUG is a publicly available code, and can be found 
at |http: / / sites.google.com/ site/runslug/, 
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Fig. 9. — i?-band, FUV, and ionizing photon luminosities vs. time for galaxies with constant SFRs of 1, 10 ^, and 10 ^ Mq yr ^ as 
indicated. i?-band and FUV luminosities are in units of erg s~^ Hz~^. We compare a fully sampled realization from SB99 (solid black lines) 
with 100, 500, and 1000 reahzations from SLUG for SFRs of 1, 10'^, and 10"^ Mq yr'^ respectively. The SLUG models are represented 
by their mean (black dash-dotted hne), median (colored dashed line) and 5-95 percentile range (filled color region). Our SLUG models 
were set to only output every 10 million years. Note that the y-axis in each panel has been chosen to match the SFR, but always spans 
the same logarithmic interval. 




Fig. 10. — Same as Figure|9l but this time made with unclustered star formation, and using lower SFR. Note the third panel of Figure 
[9] is the same SFR as the first panel of this figure. These figures were constructed with 100, 500, and 1000 realizations at SFRs of 10~^, 
10~^, and 10~* Mq yr~^ respectively. 
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Fig. 11. — Solid lines show the evolution of Q(H") and R band luminosity for individual simulations with clustered star formation with 
SFRs of 1, 10-\ and 10"^ M© yr'^. Dashed lines show the SB99 prediction. Note that the y-axis in each panel has been chosen to match 
the SFR, but always spans the same logarithmic interval. 
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Kroupa (2001) IGIMF 

T ' — I ' — I ' — I ' — I ' — I 1 — I ' — I •— 




Actual Med [Mq] Actual Med [Mq] 

Fig. 12. — The mass of the largest star in a cluster vs. that cluster's mass for clusters created by SLUG for a lKroupal 11200 II ) IMF (left) 
and the IGIMF (right). The black lines denote the analytic prediction of the maximum possible stellar mass in a cluster in the IGIMF 
model, the black dashed line notes the lower limit of the cluster mass function, and blue contours denote the location of SLUG models. 
Top panels show the maximum stellar mass as a function of the cluster mass drawn from the CMF, while bottom panels show the same 
relation relative to the sum of the masses of all stars actually populating the clusters. These two differ slightly- see section [3.21 



APPENDIX 
IMPLEMENTATION OF IGIMF 



The IGIMF theory is a statement that the SFR controls the upper cutoff of the CMF and that each cluster's mass 
changes the upper cutoff of the IMF in that cluster. 

We implement the IGIM F foll owinglWeidner et all j2010bD . We use the work of lPflamm-Altenburg fc Kroupal (|2008D , 
iWeidner fc Kroupal ()2005[) , and iWeidner et al.l ( 20041 ) to define the maximum cluster mass as 



Mecl,ma. = 84793 f ) ^ , (Al) 

\Mq yr V 

where (SFR) is the time-average SFR. Thus the SFR affects the upper cut off of the CMF. We determine the average 
star formation rate over a time interval defined by the user (fiducially 10^ yr). 

After a cluster mass has been drawn, we must adjust the upper cutoff of the IMF that we use to draw stars for 
that cluster. The r elatio n bet ween maximum stellar mass and cluster mass {mmax — Med) has been studied by 
IWeidner fc Kroupal (|2004[ ) and IWeidner et al.l H2010aD . Following their treatment, we solve a system of equations 
numerically for nimax as a function of Med ■ 

The first equation is simply a statement that the total cluster mass (Med) is the integral of the distribution of masses 
(^) integrated from the lowest to highest mass star in the cluster: 

Med = / m—dm. (A2) 
.L. ■ dm 
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The next constraint is derived based on the statement that there is only one star in the cluster with mass equal to 
i^max- Their choice of implementation of this statement is as foUowJ*^: 



1 = 



dN , 
— — dm 

dm 



(A3) 



where mmax,* is the maximum stellar mass possible. 

In the specific case of a .Kroupa (,2001.) IMF, these equations reduce to the following (taken from IWeidner fc Kroupal 
[200l . 
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- 0^3 



^(m?-- - m^-) + V"^"^ ^(^2-^.3 _ ^2-.3) (A5) 



where 



ao = +0.30, miow = 0.01 

ai = +1.30, mn = 0.08 

a2 = +2.30, mo = 1.00 

a-i = +2.35, m^aj;,* = 120 



(A6) 



We fit a 6th order polynomial to the numerical solution to find: 

6 

log^o mmax = ^ ai(logio MciY (A7) 

1=0 

where a=[ 1.449, -2.522, 2.055, -0.616, 0.0897, -0.00643, 0.000 182]. 

We then use this upper mass limit to modify the standard iKroupal ((2001) IMF to fill in the stars for the cluster. 
Figure [T^ demonstrates the result. One can see that we are accurately applying the cutoff to the IMF in the IGIMF. 
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